function [e] = fc_wk3(p);

global time pressure flow T

% function 3-el WK-model
% used to estimate C via 3-el WK model
% method: put flow1 in m-el model
% minimise difference between measured and computed pressure1
% p(1)=tVR
% p(2)=Z0
% p(3)=C

aantal=length(pressure);

fftQ=fft(flow);

R_WK3=p(1);
Z0_WK3=p(2);
C_WK3=p(3);

% application 3element WK model: compute pressure1 response of 3-el WK model to flow1

% step 1: transform Matlab fft flow1 into harmonics

fftQ(1)=fftQ(1)/aantal;


if rem(aantal,2)==0  
	for k=2:aantal/2-1
	    fftQ(k)=2*fftQ(k)/aantal;
	end;
else
   for k=2:ceil(aantal/2)-1
	    fftQ(k)=2*fftQ(k)/aantal;
	end;
end;


% step 2: compute input impedance 3-el WK model and harmonics P_WK3 

if rem(aantal,2)==0  
	for j=1:aantal/2-1
	    freq(j)=(j-1)/T;
   	 Z_WK3(j)=Z0_WK3+R_WK3/(1+i*2*pi*freq(j)*R_WK3*C_WK3);
	    fftP_WK3(j)=Z_WK3(j)*fftQ(j);
	end;
else
	for j=1:ceil(aantal/2)-1
	    freq(j)=(j-1)/T;
   	 Z_WK3(j)=Z0_WK3+R_WK3/(1+i*2*pi*freq(j)*R_WK3*C_WK3);
	    fftP_WK3(j)=Z_WK3(j)*fftQ(j);
	end;
end;


% step 3: compute matlab harmonics

fftP_WK3(1)=fftP_WK3(1)*aantal;

% check for even or uneven number

if rem(aantal,2)==0  
	for j=1:aantal/2-1
	    fftP_WK3(j+1)=fftP_WK3(j+1)*aantal/2;
	    fftP_WK3(aantal-j+1)=conj(fftP_WK3(j+1));
	end;
	fftP_WK3(aantal/2)=fftP_WK3(j+1)*aantal/2;
else
	for j=1:ceil(aantal/2)-1
	    fftP_WK3(j+1)=fftP_WK3(j+1)*aantal/2;
	    fftP_WK3(aantal-j+1)=conj(fftP_WK3(j+1));
	end;
end;


% step 4: compute inverse matlab transform

P_WK3=real(ifft(fftP_WK3));


if size(pressure)==size(P_WK3)
   verschil=(pressure-P_WK3).*(pressure-P_WK3);
else
   verschil=(pressure-P_WK3').*(pressure-P_WK3');
end;   

e=sum(verschil);



